//高斯消元法
//gauss_elimination.cpp

//用高斯消元法求解线性方程组

//1)基础概念
//线性方程组：
//a(1,1)x1 + a(1,2)x2 + ... + a(1,n)xn = b1
//a(2,1)x1 + a(2,2)x2 + ... + a(2,n)xn = b2
//......
//a(m,1)x1 + a(m,2)x2 + ... + a(m,n)xn = bm
//设x是解的一列向量(x1, x2, ..., xn)，b是右边一列系数向量(b1, b2, ..., bm)
//A是m*n的矩阵，第i行第j列元素是a(i,j)
//也可以把A看作一行向量(A1, A2, ..., An)
//其中每个元素Ai都是一列向量(a(1,i), a(2,i), ..., a(m,i))
//一般简写作Ax = b
//增广矩阵：
//将上面的线性方程组写成这样的形式：
//|a(1,1)  a(1,2)  a(1,3) ...  a(1,n) | b1 |
//|a(2,1)  a(2,2)  a(2,3) ...  a(2,n) | b2 |
//|a(3,1)  a(3,2)  a(3,3) ...  a(3,n) | b3 |
//......
//|a(m,1)  a(m,2)  a(m,3) ...  a(m,n) | bm |
//增广矩阵中左边部分是方程组中的m*n的矩阵A，右边部分是方程组中的系数向量b
//增广矩阵是一个m*(n+1)的矩阵，省略的x和列对应起来
//比如第1行指代a(1,1)x1 + a(1,2)x2 + a(1,3)x3 + ... + a(1,n)xn = b1
//线性无关，线性有关：
//设方程组中前三行方程为A1x = 0，A2x = 0，A3x = 0
//会出现这样的情况，比如添加一些常数的计算让A3可以由A1和A2转化而来
//计算这样的方程组时直接忽略A3这样可以由其他方程转化而来到方程
//那些无法相互转化的向量称为线性无关的，可以相互转化的向量称为线性有关的
//相互之间线性无关的方程也称有效的，可以被转化的方程也称无效的
//秩：
//将方程组中的无效方程删去，只留有效方程
//这时方程的个数，即A矩阵的行数即为方程组的秩，记作r
//在方程组中r <= n恒成立，若r = n则称方程组满秩
//高斯消元法：
//通过初等变化将增广矩阵变成类似于：
//|c(1,1)  c(1,2)  c(1,3) ...  c(1,r) ...  c(1,n) | b1  |
//|0	   c(2,2)  c(2,3) ...  c(2,r) ...  c(2,n) | b2  |
//|0	   0	   c(3,3) ...  c(3,n) ...  c(3,n) | b3  |
//|.....										  | ... |
//|0	   0	   0	  ...  c(r,r) ...  c(r,n) | br  |
//|0	   0	   0	  ......		   0	  | br+1|
//|0	   ......						   0	  | 0	|
//|.....										  | ... |
//上式中r为方程组的秩，n-r个无效的方程被消去，只留下r个有效方程
//1)无解情况
//若在从r+1到n行的矩阵中，右边存在系数不等于0，比如存在br+1 != 0
//则该方程组无解，因为无法找出一组x = (x1, x2, ..., xr, ..., xn)使得
//0*x1 + 0*x2 + ... + 0*xn = 0 = br+1 != 0
//2)有解情况
//从r+1到n行的矩阵中，右边没有不等于0的系数，不存在这样的br+1 != 0，则可以得到：
//c(1,1)*x1 + c(1,2)*x2 + ... + c(1,r)*xr = b1 - c(1,r+1)*xr+1 - ... - c(1,n)*xn
//			  c(2,2)*x2 + ... + c(2,r)*xr = b2 - c(2,r+1)*xr+1 - ... - c(2,n)*xn
//......
//								c(r,r)*xr = br - c(r,r+1)*xr+1 - ... - c(r,n)*xn
//当方程组的秩r < n，方程组有无穷解
//当方程组的秩r = n，方程组有唯一解，有效方程组经过高斯消元之后是这样的：
//|c(1,1)  0	   0	  ...  0	  | b1  |
//|0	   c(2,2)  0	  ...  0	  | b2  |
//|0	   0	   c(3,3) ...  0	  | b3  |
//|.....							  | ... |
//|0	   0	   0	  ...  c(r,r) | br  |
//其中的r = n
//r = n，即方程组的秩等于未知量的个数的方程组有唯一解，通过高斯消元法可以求出
//具有无穷解的方程组的求解需要确定一组基础解系，比较麻烦
//本文中只求解具有唯一解的线性方程组
//高斯若尔当消元法：
//Gauss-Jordan消元法是高斯消元法的另一个版本，该算法最终可以得到一个这样的矩阵：
//|c(1,1)  0	   0	  ...  0	  ...  c(1,n) | b1  |
//|0	   c(2,2)  0	  ...  0	  ...  c(2,n) | b2  |
//|0	   0	   c(3,3) ...  0	  ...  c(3,n) | b3  |
//|.....										  | ... |
//|0	   0	   0	  ...  c(r,r) ...  c(r,n) | br  |
//|0	   0	   0	  ......		   0	  | br+1|
//|0	   ......						   0	  | 0	|
//|.....										  | ... |
//从而可以直接得到x1到xr每个未知量的解，没有高斯消元法中最后的向上回溯迭代的过程
//Gauss-Jordan消元法的时间复杂度较高
//
//2)具体实现
//在代码的实现中，将方程组通过消元和代入两个过程进行处理，最终得到方程组的解
//消元和代入两个操作是对数学行为的代码化
//消元和代入在线性代数课程中有详细的讲解，还有许多关于行列式和矩阵的各种性质与公式
//由于篇幅限制本文不再详述
//更多内容见“数学复习全书(2013年李永乐李正元考研数学 数学一)”，作者“李永乐”“李正元”
//
//本文引用了“C#算法应用之高斯消元法实现”，作者“佚名”
//“高斯消元法(Gauss Elimination)分析&题解&模板——czyuan原创”，作者“Keep Moving 永不止步...”

#include "general_head.h"
#include "matrix.h"
double select_pivot(matrix e, int n, int p, int& index);
void elimination(matrix& e, int n);
void backward_substitution(matrix e, int n, double *x);
void forward_substitution(matrix e, int n, double *x);

void gauss_elimination(matrix e, int n, double *x)
{//满秩的增广矩阵e是n行n+1列的，行下标从0到n-1，列下标从0到n
 //其中从0到n-1列是线性方程组中的矩阵A，最后一列是线性方程组中的系数向量b
 //返回该线性方程组的唯一解，存储于数组x中，数组x长度为MAX
	elimination(e, n);
	backward_substitution(e, n, x);
}
void elimination(matrix& e, int n)
{//增广矩阵e是n行n+1列的，最后一列是系数向量b
 //通过消元将增广矩阵的左边部分化简为一个上三角矩阵
 //满秩方程组中不会出现某一列系数全为0的情况
 //假如设增广矩阵中i列全为0，则未知量xi可以取任意值从而该方程组存在无穷解
	for(int k = 0; k < n - 1; ++ k){
		//k是n阶矩阵中的主对角线，n阶矩阵的消元是沿着主对角线进行的
		int index;
		double pivot = select_pivot(e, n, k, index);
		//若index并不是起始行k，则交换index行和k行
		if(pivot != e.m_m[k][k])
			//交换从k列开始，保证之前从0到k-1列元素都已经是0
			for(int i = k; i <= n; ++ i)
				swap(e.m_m[index][i], e.m_m[k][i]);
		//用k行k列的系数消除从k+1行到最后一行的k列系数
		//保证之后所有行中，从0到k列的系数全为0
		for(int i = k + 1; i < n; ++ i){
			//从k+1到最后一行
			//tmp是i行k列元素除以k行k列元素的商
			double tmp = e.m_m[k][k] / e.m_m[i][k];
			for(int j = k; j <= n; ++ j)
				//使i行k列元素为0，且该行之后的所有元素也满足等式
				e.m_m[i][j] = e.m_m[i][j] * tmp - e.m_m[k][j];
		}
	}
}
double select_pivot(matrix e, int n, int k, int& index)
{//增广矩阵e是n行n+1列的，对k行k列进行选择主元的操作
 //返回替换行的行标号index和主元的绝对值
	double pivot(fabs(e.m_m[k][k]));
	index = k;
	//找出从k行开始到最后一行中，k列上元素的绝对值最大的一行index
	//保证从k行到最后一行，从0到k-1列的元素都是0，系数都已被消去
	for(int i = k + 1; i < n; ++ i)
		//考察的是元素的绝对值
		if(pivot < fabs(e.m_m[i][k])){
			pivot = fabs(e.m_m[i][k]);
			index = i;
		}
	return(pivot);	
}
void backward_substitution(matrix e, int n, double *x)
{//秩为n的方程组对应的增广矩阵e是n行n+1列的，最后一列是系数向量b
 //增广矩阵e中左边的n阶矩阵是一个这样的上三角矩阵：
 //| a(1,1)  a(1,2)  a(1,3)  a(1,4)  ...	a(1,n) | b1 |
 //| 0	     a(2,2)  a(2,3)	 a(2,4)	 ...	a(2,n) | b2 |
 //| 0		 0		 a(3,3)  a(3,4)	 ...	a(3,n) | b3 |
 //......
 //| 0		 0		 0		 0	     ...	a(n,n) | bn |
 //由最后一行可以直接通过a(n,n)*xn = bn计算出xn = bn / a(n,n)
 //然后向上依次代入每个方程即可求出x的唯一解x = (x1, x2, x3, ..., xn)
 //这个求解过程称为逆向代入
	for(int i = n - 1; i >= 0; -- i){
		//满秩方程组有唯一解的情况中
		//方程组中的矩阵A的行数与列数与方程组的秩，三个数相等
		x[i] = e.m_m[i][n] / e.m_m[i][i];
		for(int j = i; j >= 0; -- j){
			//从i行逆向向上代入每一行
			//j行n列，即j行的b系数减去xi乘以j行i列的系数
			e.m_m[j][n] -= e.m_m[j][i] * x[i];
			//j行i列的系数置0
			e.m_m[j][i] = 0;
		}
	}
}
void forward_substitution(matrix e, int n, double *x)
{//秩为n的方程组对应的增广矩阵e是n行n+1列的，最后一列是系数向量b
 //增广矩阵e中左边的n阶矩阵是一个这样的下三角矩阵：
 //| a(1,1)	 0		 0		 0	   	 ...	0	   | b1 |	
 //| a(2,1)  a(2,2)  0		 0		 ...	0	   | b2 |
 //| a(3,1)  a(3,2)  a(3,3)  0       ...	0	   | b3 |
 //......
 //| a(n,1)  a(n,2)  a(n,3)  a(n,4)  ...	a(n,n) | bn |
 //由第1行可以直接通过a(1,1)*x1 = b1计算出x1 = b1 / a(1,1)
 //然后向下依次代入每个方程即可求出x的唯一解x = (x1, x2, x3, ..., xn)
 //这个过程称为正向代入
	for(int i = 0; i < n; ++ i){
		x[i] = e.m_m[i][n] / e.m_m[i][i];
		for(int j = i; j < n; ++ j){
			e.m_m[j][n] -= e.m_m[j][i] * x[i];
			e.m_m[j][i] = 0;
		}
	}
}
